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ABSTRACT 

A  transient  three-dimensional  full- field  analysis  of  a  three-point-bend  fracture  specimen  dy¬ 
namically  loaded  into  the  fully  yielded  plastic  state  has  been  carried  out.  The  specimen  is  rapidly 
loaded  by  means  of  a  concentrated  transverse  force  applied  at  mid-span  on  the  uncracked  surface  of 
the  specimen.  It  is  assumed  that  the  material  is  ductile  and  fracture  initiation  occurs  after  substan¬ 
tial  plastic  deformation  has  developed  in  the  uncracked  ligament.  Furthermore  it  is  supposed  that 
the  crack  tip  conditions  are  such  that  the  J-integral  may  be  adopted  as  a  characterizing  parameter. 
We  derive  an  expression  for  the  local  energy  flux  appropriate  to  a  three-dimensional  crack  front. 
Based  on  this  fundamental  crack  tip  flux  integral,  a  domain  integral  representation  for  J  which 
is  naturally  suited  for  the  finite  element  method  is  obtained.  Using  the  domain  integral,  accurate 
point  wise/local  values  of  J  along  a  three-dimensional  crack  front  can  be  extracted  from  numerically 
determined  field  quantities.  The  effect  of  transient  loading,  geometry  and  plastic  deformation  on 
the  variation  of  J  along  the  crack  front  is  examined.  A  purpose  of  the  present  study  is  to  determine 
conditions  under  which  the  value  of  J  at  initiation  may  be  inferred  from  quantities  that  are  directly 
measurable  in  a  dynamic  fracture  experiment.  In  an  earlier  paper  a  transition  time  was  introduced 
to  provide  a  practical  bound  on  the  time  range  over  which  conditions  of  J-dominance  prevailed 
and  the  deep-crack  formula  was  applicable  under  transient  loadings.  The  transition  time  concept  is 
further  examined  in  this  paper.  With  the  full-field  solutions  in  hand,  we  comment  on  the  suitability 
of  proposed  surface  locations  for  measurements  of  moment  and  rotation,  the  input  measurements 
for  application  of  the  deep- crack  formula  for  J.  Implications  for  fracture  testing  of  tough  materials 
at  relatively  high  loading  rates  are  discussed. 


*currently  at  Department  of  Mechanical  Engineering,  M.I.T.,  Cambridge,  MA  02139 
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1.  INTRODUCTION 


Several  experimental  procedures  have  been  recently  proposed  for  the  measurement  of  dynamic 
fracture  toughness.  In  an  earlier  study  [1],  we  carried  out  an  analysis  of  a  three-point-bend  frac¬ 
ture  specimen  under  dynamic  loading,  modeling  a  specimen  that  Joyce  and  Hackett  have  used  to 
determine  the  fracture  toughness  of  ductile  structural  steel  [2].  In  [1],  a  full-field  plane  strain  finite 
element  analysis  is  described.  The  crack  tip  values  of  J  were  computed  by  means  of  a  domain 
integral  expression  and  the  results  were  compared  with  the  estimate  of  J  determined  from  the 
deep-crack  formula  proposed  by  Rice,  Paris,  Merkle  [3]  (with  parameters  appropriately  replaced 
for  the  dynamic  analysis).  Furthermore,  a  notion  of  transition  time  was  introduced  to  provide  a 
practical  bound  on  the  time  range  over  which  conditions  of  J-dominance  prevailed  near  the  crack 
front  and  the  deep-crack  formula  was  applicable  under  transient  loadings. 

In  order  to  completely  understand  the  response  of  a  fracture  test  specimen,  a  three-dimensional 
analysis  is  required.  This  is  especially  important  if  there  is  substantial  plastic  flow  in  the  ligament 
of  the  specimen.  Several  analytical  investigations  based  on  three-dimensional  models  have  been 
carried  out  for  the  case  of  quasi-static  loading  [4-6].  In  these  studies,  the  J-integral  [7]  is  viewed  as 
a  paramenter  characterizing  the  crack  tip  field  and  the  determination  of  J  as  a  function  of  position 
along  crack  front  was  a  principal  objective.  The  studies  revealed  that  the  variation  of  J  along  crack 
front  is  strongly  dependent  on  the  amount  of  plastic  deformation  in  the  uncracked  ligament  of  the 
specimen. 

A  full-field  three-dimensional  finite  element  analysis  of  a  dynamically  loaded  three-point-bend 
fracture  specimen  is  presented  in  this  paper.  First  a  general  expression  for  the  crack  tip  J-integral 
based  on  the  fundamental  energy  flux  integral  is  given  for  a  three-dimensional  crack  front  in  a 
deformable  body  subject  to  transient  loadings.  Then,  a  domain  integral  form  of  J  is  derived  from 
the  crack  tip  integral,  and  an  expression  is  given  for  pointwise  value  of  J  along  a  crack  front. 
Using  the  domain  integral  representation,  accurate  J  values  can  be  extracted  from  finite  element 
solutions.  The  domain  integral  formulation  corresponds  to  the  method  of  virtual  crack  extension 
[8-11]  and  is  particularly  suited  for  evaluation  of  J  in  three-dimensional  crack  problems.  For  the 
numerical  simulation,  a  model  of  the  three-point-bend  fracture  specimen  is  constructed  with  three- 
dimensional  finite  elements.  The  transient  finite  element  analysis  is  carried  out  for  the  time  interval 
of  interest  using  an  explicit  time  integration  scheme.  Pertinent  quantities,  including  J  as  a  function 
of  position  along  the  crack  front,  are  computed  and  recorded  during  the  calculations.  From  the 
results,  interpretations  are  made  on  the  J  variation  along  the  crack  front  and  its  dependence  on 
loading  rate  and  elapsed  time.  Strong  variation  of  field  variables  in  the  through-the-thickness 
direction  is  observed  and  these  three-dimensional  effects  are  assessed. 
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In  view  of  the  experimental  procedure,  the  accuracy  of  J  estimated  from  the  deep-crack  formula 
is  evaluated  by  comparing  estimated  values  with  precise  values  extracted  from  the  full-field  solution 
using  the  domain  integral  expression.  In  addition,  we  will  comment  on  the  suitability  of  proposed 
surface  locations  for  measurements  of  moment  and  rotation,  these  being  the  input  measurements  for 
appplication  of  the  deep-crack  J  formula.  The  deep-crack  expression  is  essentially  a  two-dimensional 
formula  and  we  will  discuss  its  applicability  in  the  three-dimensional  fracture  problems  under  both 
quasi-static  and  dynamic  loading  conditions.  The  present  three-dimensional  analysis  clarifies  these 
and  other  geometrical  effects  associated  with  a  standard  three-point-bend  fracture  specimen. 
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2.  DYNAMIC  /-INTEGRAL 


2.1  Energy  integral  for  three-dimensional  crack  front 


In  the  analysis  of  dynamic  crack  growth  in  an  elastic  solid  under  two-dimensional  conditions, 
the  dynamic  energy  release  rate  J  is  defined  as  the  energy  released  from  the  body  per  unit  crack 
advance.  More  precisely,  it  is  the  limiting  value  of  the  energy  flux  across  a  path  T  which  surrounds 
the  crack  tip,  divided  by  the  instantaneous  crack  tip  speed,  as  the  path  is  shrunk  onto  the  crack 
tip.  A  mathematical  expression  for  the  energy  release  rate  in  terms  of  the  crack  tip  fields  for 
two-dimensional  problems  with  crack  extending  in  the  X\ -direction  [12]  is 


J  —  lim 
r— o 


L 


dxt ' 

( U  +  T)n  1  - 


dT 


(2.1) 


Here  <7,^-  and  u,  are  the  cartesian  components  of  stress  and  displacement,  and  n;  are  the  components 
of  a  unit  vector  normal  to  T  and  pointing  away  from  the  crack  tip.  The  quantities  U  and  T  are  the 
stress  work  and  the  kinetic  energy  densities,  respectively,  and  are  defined  as 


U  = 


d2Ui 


dt  'dxj 


dt' 


1  dui  dui 
~  2P  dt  dt 


(2.2) 


To  extend  this  idea  to  three-dimensional  fields,  suppose  that  the  crack  edge  is  defined  in 
rectangular  coordinates  by  £i(<s,t)  at  any  time  £,  where  s  is  the  arclength  along  the  crack  edge 
measured  from  some  arbitrary  point.  The  speed  of  the  crack  edge  at  any  s  is  v(s)  =  and 

the  direction  of  crack  advance  is  locally  —  £{/v  (see  Fig.  la). 

Consider  a  tubular  surface  St  enclosing  the  crack  edge  and  moving  with  it.  The  tubular  surface 
is  formed  as  follows.  In  any  plane  locally  perpendicular  to  the  crack  edge,  specify  a  small  crack  tip 
contour  T  that  begins  on  one  face  of  the  crack,  encircles  the  crack  edge,  and  ends  on  the  opposite 
crack  face.  The  tubular  surface  is  completely  specified  by  the  condition  that  its  intersection  with 
every  plane  locally  perpendicular  to  the  crack  edge  is  the  same  curve  T  at  any  time  t  (see  Fig.  lc). 

For  a  crack  advancing  under  three-dimensional  conditions,  at  any  instant  of  time  the  energy 
flux  through  the  tubular  surface  St  is 


(  U  T  T  &ijTlj'U'i 


dSt 


(2.3) 


where  n{  are  the  components  of  the  unit  normal  vector  to  St  directed  away  from  the  crack  edge  and 
other  quantities  are  as  defined  above.  At  any  point  on  5*,  the  particle  velocity  U{  may  be  separated 
into  two  contributions,  one  arising  from  the  fact  that  U{  is  changing  at  the  fixed  point  on  St  and  a 
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convected  contribution  arising  from  the  fact  that  this  point  is  moving  through  material  where  it; 
has  a  spatial  gradient.  Thus, 


Ui  = 


diii 

~dt 


_  du±: 

st  d*i 


Therefore,  if  $a  <  s  <  s*,  is  the  range  of  arclength  along  the  crack  edge,  then 


F(St)  =  I*  v(a)  jf 


du{ 

(■ U  +  T)vknk  -  —  uk 


dT  ds 


(2.4) 


(2.5) 


where  dT  is  that  part  of  the  tube  between  s  and  s  +  ds.  The  asymptotic  relation,  uf  ~  -(du,i/dxj)£j 
as  r  — »  0+,  has  been  invoked  in  (2.5)  [12].  The  quantity  v(s)  is  the  local  crack  advance  per  unit 
time,  so  that  the  integral  multiplying  it  is  the  energy  flow  per  emit  crack  advance  through  dT .  The 
local  energy  release  rate  is 


J(s) 


vk  lim 
r  * 


im  / 

-*°  Jr 


{U  T)nk  &ijHj 


dxtj 

dxk 


dr 


(2.6) 


The  direction  of  crack  advance  vk  depends  only  on  s  and  is  taken  out  of  the  integration.  Under 
quasi-static  conditions  and  where  U  is  taken  to  be  the  strain  energy  function,  the  integral  to  the 
right  of  vk  is  the  so-called  Jk  conservation  integral  [13,14].  A  unified  treatment  of  crack  tip  integrals 
and  conservation/dissipation  integrals  as  direct  consequences  of  appropriate  balance  laws  can  be 
found  in  [15]. 

Suppose  that  over  an  arbitrarily  small  increment  of  time,  the  crack  front  at  s  advances  by 
A(s)  in  the  normal  direction  (in  the  plane  of  the  crack)  within  the  segment  sa  <  s  <  With 
respect  to  the  reference  coordinate  system,  the  crack  growth  increment  is  A (s)vk($)  which  is  more 
conveniently  denoted  by  lk(s)  as  depicted  in  Fig  lb.  Associated  with  this  increment  of  growth  of 
crack  front,  the  total  energy  released  by  the  body  is 

fSb 

J  =  /  J(s)X(s)ds=  lim 
/  r-o 


(U  T  T)nk 


dui 

dxk 


h  dS 


(2.7) 


where  dS  =  dVds.  The  energy  released  through  different  segments  of  the  crack  front  is  obtained 
from  (2.7)  by  appropriate  choice  of  sa  and  s^. 

For  quasi-static  monotonic  loading  conditions  and  proportional  stress  history  at  each  material 
point,  the  integral  (2.1)  is  path  independent  for  a  nonlinear  elastic  or  elastic-plastic  material,  and 
the  value  of  J  does  not  depend  on  the  limiting  process  of  shrinking  the  contour  V  onto  the  crack 
tip.  In  this  case,  the  integral  is  precisely  Rice’s  /-integral  [7]  and  its  value  is  the  amplitude  of 
near  crack  tip  field.  Under  dynamic  loading  conditions,  inertial  effects  preclude  the  definition  of  an 
equivalent  path  independent  integral  and  the  value  of  the  near  tip  field  intensity  must  be  expressed 
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in  terms  of  the  crack  tip  limit  in  (2.1),  although  the  shape  of  the  contour  as  shrunk  onto  the  crack 
tip  is  arbitrary.  For  a  3-D  crack  front,  the  local/pointwise  value  of  J  is  only  defined  in  terms  of  the 
limit  in  (2.6). 

2.2  Domain  integral  representation  for  J 

It  is  apparent  that  from  a  discrete  computational  point  of  view,  the  expression  (2.6)  is  not 
suitable  for  evaluating  values  of  J(s )  nor  is  (2.7)  useful  for  evaluating  values  of  J .  As  in  the  two- 
dimensional  case,  accurate  evaluations  of  field  near  the  crack  front  is  difficult  and  alternate  forms 
for  J  that  are  better  suited  for  numerical  calculations  are  discussed  here.  The  finite  domain  form 
for  J  follows  directly  from  the  fundamental  crack  tip  energy  integral  (2.7)  by  the  application  of  the 
divergence  theorem  over  part  of  the  cracked  solid.  The  present  derivation  follows  the  formulation 
in  [6,11]  where  it  is  presented  in  greater  detail. 

Let  S0  be  an  arbitrary  outer  surface  (including  the  end-caps)  which  surrounds  the  inner  tubular 
surface  St  as  depicted  in  Fig.  lc.  The  connecting  upper  and  lower  crack  faces  are  S+  and  S- .  Thus 
S  (S  =  S0  +  St  +  S+  +  S~)  is  a  closed  surface.  To  write  (2.7)  as  an  integral  over  a  closed  and  an 
integral  over  physical  boundaries  (i.e.  <S_j_  -f  <S_)  we  introduce  a  weighting  function  q k.  The  vector 
qk  is  a  smooth  function  of  x\,  and  x$  subject  to  these  restrictions.  It  is  equal  to  /^(s)  on 
vanishes  on  S0  and  varies  in  an  appropriate  manner  (to  be  specified)  along  the  crack  faces.  With 
qk  so  defined,  the  integral  in  (2.7)  can  be  written  as  an  integral  over  the  closed  surface  S  and  an 
integral  over  the  physical  boundaries  <S+  +  <S_, 


~(U  +  T)Sjk 


mjqk  dS 


-  lim  f 
r^°  Js++S- 


dui 


afk-(u  +  T)Sli 


mjqk  dS 


(2.8) 


Here  raj  is  the  unit  normal  on  S  that  points  away  from  the  enclosed  volume.  On  the  surface  St , 
?7i j  is  the  negative  of  nj  which  has  been  defined  earlier.  On  each  of  the  crack  faces,  qk  is  prescribed 
to  be  orthogonal  to  mk  (i.e.  qkmk  =  0  on  and  S- ).  In  addition  if  crack  faces  are  traction  free, 
the  integral  over  S+  +  S„  in  (2.8)  vanishes  identically. 

We  apply  the  divergence  theorem  to  the  closed  integral  in  (2.8)  and  invoke  the  balance  of  linear 
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momentum  to  get  the  following  finite  domain  representation 


(2.9) 


In  the  above  equation,  crack  faces  are  assumed  to  be  traction  free.  If  crack  face  traction  cannot 
be  neglected,  the  second  integral  in  (2.8)  must  be  appended  to  the  integral  in  (2.9).  The  inclusion 
of  traction  on  crack  faces  and  body  force  like  terms  are  discussed  in  greater  detail  in  [6].  In  (2.9), 
the  domain  V  has  been  taken  to  be  the  total  volume  enclosed  by  the  arbitrary  outer  surface  <S0, 
including  the  crack  tip  region.  Any  difference  between  the  integration  over  the  volume  V  and  the 
volume  enclosed  within  S  vanishes  as  St  is  shrunk  onto  the  crack  front.  It  may  be  noted  that  the 
last  two  terms  in  (2.9)  cancel  identically  in  the  case  of  an  elastic  solid. 

We  emphasize  that  any  smooth  (C°)  function  which  is  appropriately  defined  within  the  enclosed 
volume  (i.e.  the  function  assume  the  precribed  values  on  surfaces  <St,  SQ  and  crack  faces  <S+,  S-) 
could  serve  as  the  weighting  function  qk-  The  vector  function  qk{xux2‘>xz)  maY  be  interpreted  as 
the  virtual  translation  of  material  point  (aq,  x*i ,  xs)  in  V  due  to  virtual  extension  of  the  crack  front 
by  /^(s)  [6,8-11].  Since  SQ  is  an  arbitrary  outer  surface,  the  value  of  J  does  not  depend  on  the  size 
or  shape  of  the  arbitrary  volume  V  enclosed  by  Sa.  This  property  of  the  domain  representation 
provides  a  useful  check  on  the  consistency  and  quality  of  the  numerical  solution. 

An  approximation  of  the  value  of  J(s)  is  obtained  by  assuming  that  it  is  constant  in  an  interval 
sa  <  s  <  Sfc.  The  crack  front  is  then  perturbed  in  this  interval,  say  an  amount  A (s)  in  the  plane  of 
the  crack,  where  A(^)  —  0  outside  of  the  interval.  This  perturbation  of  the  crack  front  results  in  a 
total  energy  release  rate  J  according  to  (2.9).  Then  J(s)  is  given  by 

J(s)  -1  j  l  a  A (s)ds  (2.10) 

J  3b 

A  more  precise  scheme  for  calculating  J(s)  which  is  also  consistent  with  the  finite  element  formu¬ 
lation  is  discussed  in  Section  4.2. 

2.3  Specialization  to  a  stationary  crack  front 

For  the  case  of  stationary  crack  tip  in  a  dynamically  loaded  body,  the  kinetic  energy  is  bounded 
at  the  crack  tip  and  consequently  it  makes  no  contribution  to  (2.1).  The  energy  release  rate  for 


-L 


dui  dqk  ,TT  ,  dqk  ,  . 
aij  —  —  ~  (U  +  T)~  +  P 


dxk  dxj 


dxk 


d2Ui  dui  du{  d2U{ 
dt2  dxk 


+ 


dt  dxkdt 

d2U{ 

dxjdxk 


Qk 


dU  \ 
dxk)qk 


dV 
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this  case  is 


J  =  fisi /r  (r,n'  -  (2-u) 

Similarly  the  expression  for  the  local  energy  release  rate  at  the  point  s  on  the  three-dimensional 
crack  is 

J{s)  =  lim  f  (unk  -  <Tijnj~-\  ukdT  (2.12) 

r^0JT{s)  V  dxkJ 

The  specialization  of  (2.9)  to  a  stationary  crack  front  in  a  dynamically  loaded  body  is 


J  = 


duj  dqk 
dxk  dxj 


d2Ui  du{ 

+  P~dt2  dxk 


dV 


(2.13) 


It  may  be  noted  that  (2.13)  can  be  derived  directly  from  (2.12)  using  the  procedure  which  led  to  the 
integral  in  (2.9).  A  method,  based  on  the  expression  (2.13),  for  extracting  accurate  values  of  the 
intensity  of  near  tip  stresses  and  deformation  from  the  fields  determined  from  numerical  simulation 
will  be  presented  in  Section  4. 

To  arrive  at  the  expression  (2.13),  U  is  taken  to  be  a  potential  function  for  stress,  i.e., 


dU 


aij  8e 


(2.14) 


Here  eij  are  the  cartesian  components  of  strain.  In  the  case  of  a  solid  whose  mechanical  response  is 
described  by  a  stress-strain  relationship  of  incremental  type,  the  stress  work  density  is  dependent 
on  strain  history  and  (2.14)  generally  does  not  apply.  However,  when  certain  conditions  are  met 
any  deviation  from  the  equality  is  very  small  and  can  be  neglected.  The  conditions  require  the 
load  applied  to  the  crack  tip  region  to  be  monotonically  increasing  and  stress  histories  for  material 
particles  to  be  nearly  proportional  [16]. 
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3.  J  ESTIMATES  FOR.  TRANSIENT  PROBLEMS 
3.1  Transition  time 


The  dynamically  loaded  three-point-bend  specimen  has  been  proposed  for  the  measurement  of 
dynamic  fracture  toughness  [2].  A  schematic  of  the  specimen  is  shown  in  Fig.  2a.  The  specimen 
has  a  through-thickness  planar  crack  of  length  a  (the  crack  front  being  parallel  to  the  a3-axis)  and 
is  supported  by  rollers  seperated  by  a  distance  of  2 L.  A  dynamic  load  P  is  applied  on  the  surface 
Ei  =  0  at  the  center  span.  In  [1],  we  introduced  a  transition  time  tx  to  provide  an  estimate  of  the 
time  beyond  application  of  the  loading  at  which  a  /-dominated  field  is  established  in  the  crack  tip 
region  and  a  deep- crack  /-formula  is  applicable.  We  assumed  that  an  estimate  of  the  transition 
time  can  be  determined  from  the  time  history  of  the  relative  magnitudes  of  the  total  kinetic  energy 
and  the  total  deformation  energy  (or  stress  work)  of  the  specimen.  In  particular,  the  transition 
time  is  defined  as  the  time  beyond  which  the  kinetic  energy  of  the  specimen  is  less  than  the  energy 
of  deformation. 


Direct  measurement  of  the  total  kinetic  and  deformation  energies  in  a  laboratory  specimen 
is  not  possible.  However,  these  quantities  can  be  approximated  by  considering  simple  models 
consistent  with  actual  boundary  measurements.  In  [1]  an  estimate  of  the  kinetic  energy  was  obtained 
by  means  of  a  model  based  on  elastic  Bernoulli- Euler  beam  theory  with  an  assumption  that  the 
kinetic  energy  at  the  early  stage  is  dominated  by  elastic  structural  response.  To  approximate  the 
deformation  energy,  a  quasi-static  elastic  2-D  three-point-bend  model  was  considered.  For  this 
model,  the  relationship  between  the  total  strain  energy  of  the  specimen  and  the  displacement  at 
the  load  point  is  known.  Using  these  models,  the  ratio  of  kinetic  energy  to  deformation  energy 
at  time  t  is  given  in  terms  of  the  deflection  A (£)  and  the  velocity  A(t)  at  the  load-point;  these 
quantities  are  of  particular  interest  since  they  can  be  measured  in  a  dynamic  fracture  experiment. 
Specifically, 


#)  -■* 
W  /  model 


h  A(oy 

c0  A (t)J 


(3.1) 


Here  H  is  the  width  of  the  specimen,  c0  is  the  sound  speed  in  the  specimen  (i.e.,  longitudinal  bar 
wave  speed)  and  the  (time-independent)  dimensionless  shape  factor  S  depends  on  the  dominant 
mode  shape  and  elastic  compliance  Cs  of  the  specimen.  For  a  bend  specimen  of  thickness  B  and 
span  2 L  between  supports  as  shown  in  Fig.  2a,  S  is  given  by 


S  = 


LB  E  Cs 
H 


(3.2) 


where  E  is  the  Young’s  modulus.  The  values  of  S  for  a  standard  ASTM  three-point-bend  specimen 
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is  given  in  [1].  In  the  above  model  analysis  the  effect  of  discrete  stress  waves  was  ignored.  Obviously, 
the  values  of  estimated  energies  are  not  valid  until  after  elastic  waves  make  several  passes  within 
the  specimen.  The  limitation  is  of  no  consequence  since  such  wave  effects  are  unimportant  at  the 
transition  time.  The  model  analysis  also  assumed  that  effect  of  plastic  deformation  is  not  significant 
prior  to  the  transition  time. 


For  the  purpose  of  writing  an  explicit  expression  for  t.T  we  introduce  a  dimensionless  displace¬ 
ment  coefficient  D  defined  by 


D  = 


t  A (t) 


(3.3) 


tr 


For  example,  if  the  time  variation  of  the  displacement  can  be  represented  by  A  =  f3  fl  then  D  —  7. 
With  tT  defined  as  the  time  when  the  ratio  K/W  is  equal  to  unity,  we  combine  (3.1)  and  (3.3)  to 
[1] 

tT  =  D  S  —  (3.4) 

Cn 


To  evaluate  the  accuracy  of  the  estimate  of  the  transition  time  based  on  (3.1)  the  finite  element 
results  for  the  dynamically  loaded  bar  is  employed.  The  details  of  the  finite  element  model  and 
calculations  are  given  in  Section  5.  With  A (t)  and  A (t)  obtained  from  the  results  of  the  3-D  finite 
element  calculations,  and  values  of  5,  II  and  c0  appropriate  to  the  specimen,  we  calculated  the  time 
variation  of  the  ratio  K/W  according  to  (3.1).  The  ratio  of  the  actual  total  kinetic  and  deformation 
energies  from  the  finite  element  model  (obtained  by  summing  the  respective  energies  over  all  the 
elements)  is  also  calculated.  Both  ratios  are  plotted  in  Fig.  3.  The  close  agreement  between  the 
curves  is  evidence  of  the  accuracy  of  the  estimate  of  relative  energy  levels  based  on  the  elementary 
model  analysis  from  which  (3.1)  which  was  derived.  During  the  early  stages  of  the  transient,  the 
kinetic  energy  dominates.  With  progress  of  time  deformation  energy  (stress  work)  dominates.  The 
ratio  K/W  computed  by  (3.1)  decreases  to  unity,  at  time  tTCi/H  ~  28  (or  tTC0/H  ~  24)  in 
dimensionless  units.  Here  C\  is  the  longitudinal  wave  speed  in  an  unbounded  medium  and  H / C\  is 
the  time  it  takes  the  wave  to  travel  the  specimen’s  width.  It  is  worth  noting  that  the  transition  time 
determined  from  a  2-D  analysis  is  txCx/H  a  27  [1].  We  point  out  that  nearly  identical  estimates  of 
the  transition  time  are  obtained  upon  replacing  A(£)/A(t)  by  #l(0/#l(0  in  (3.1)  where  #z,(£)  is 
the  specimen  rotation  about  a  hinge  axis  at  the  ligament.  The  precise  definition  of  Oi(t)  is  given 
in  the  next  section. 


In  summary,  we  observe  that  the  response  of  specimen  can  be  conveniently  characterized 
by  a  short-time  response  dominated  by  discrete  waves  and  a  long-time  response  dominated  by 
deformation  energy.  At  the  transition  time  structual  inertial  effects  are  important.  Beyond  the 
transition  time,  inertial  effects  diminishes  relative  to  the  overall  energy  absorbed  by  the  body. 
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3.2  The  deep-crack  estimation  procedure 

On  the  basis  of  transient  2-D  finite  element  analyses  a  formula  for  the  computation  of  dynamic 
J  from  measurable  quantities  was  proposed  in  [1].  The  formula  is  a  modified  version  of  a  deep-crack 
formula  for  calculating  the  value  of  J  under  essentially  equilibrium  conditions  [3].  Under  high-rate 
of  loading,  the  inertial  resistance  of  the  specimen  screens  the  crack  tip  region  from  the  applied 
loads.  To  minimize  this  effect,  the  variables  in  the  quasi-static  formula  are  replaced  by  equivalent 
variables  which  characterize  the  near  crack  region  of  the  body.  Thus  the  moment  is  taken  to  be  the 
net  moment  Ml  carried  by  the  ligament,  and  the  corresponding  rotation  is  replaced  by  the  crack 
mouth  opening  displacement  divided  by  the  distance  between  the  crack  mouth  and  the  hinge  axis 
on  the  ligament  and  is  denoted  by  Ol ■  With  these  changes,  the  formula 

Jdc{t)  =  ~j  ML(t')deL(t')  (3.5) 

is  proposed  for  estimating  the  value  of  J  in  a  dynamically  loaded  three-point-bend  ductile  fracture 
specimen.  In  the  above  expression,  b  is  the  ligament  length  and  B  is  the  specimen  thickness.  At 
the  level  of  beam  approximation,  the  integral  of  the  moment  and  angle  in  (3.5)  represents  the  work 
done  on  the  ligament.  We  employ  the  subscript  dc  on  J  to  distinguish  estimates  based  on  the 
assumption  of  deep-crack  from  precise  values  based  on  (2.13). 

3.3  Inferred  rotation  and  moment  at  the  ligament 

To  calculate  Jdc{t)  according  to  (3.5),  the  values  of  Oi(t)  and  Mi{t)  must  be  known.  In  an 
experiment,  the  angle  9 l  can  be  determined  from  the  measurement  of  the  crack  mouth  opening 
displacement  S  and  an  estimate  of  the  hinge  axis  location  necessary  to  convert  this  displacement 
to  a  rotation.  The  opening  displacement  at  the  crack  mouth  is  essentially  constant  through  the 
thickness.  Thus  6  is  a  well-defined  measurable  quantity.  The  hinge  axis  is  an  effective  line  on 
the  ligament  plane  (and  parallel  to  the  x3-axis)  where  the  axial  strain  (and  stress)  vanishes.  The 
determination  of  the  hinge  axis  from  the  full-field  solution  will  be  taken  up  in  Section  4.3.  Under 
fully  yielded  conditions,  the  hinge  axis  can  be  rather  accurately  estimated  using  the  Green  and 
Hundy  plane  strain  slip-line  solution  [17]  for  a  rigid-perfectly  plastic  material.  The  rotation  about 
the  ligament  based  on  the  latter  hinge  axis  is  denoted  by  6*L, 

The  direct  measurement  of  the  net  moment  Ml  carried  by  the  ligament  is  extremely  difficult, 
due  in  part  to  severe  plastic  deformation  between  the  load-point  and  the  crack  tip.  Furthermore 
the  three-dimensional  nature  of  the  deformation  fields  near  the  crack  tip  and  in  the  vicinity  of  the 
ligament  preclude  the  use  of  surface  measurements  in  these  regions  for  the  determination  of  the  net 
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moment  carried  by  the  uncracked  ligament.  Thus  it  appears  the  instantaneous  values  of  ML  must 
be  inferred  from  measurements  of  the  moment  and  the  shear  force  at  some  distance  away  from  the 
cracked  section  and  from  the  roller  support.  At  a  location  remote  from  the  cracked  section  and 
the  roller  support,  (e.g.  the  section  mid-way  between  the  crack  plane  and  the  roller  support),  the 
section  should  remain  elastic  and  the  fields  should  not  vary  through  the  thickness.  Indeed  results 
from  3-D  finite  element  calculations  which  are  discussed  in  Section  5  support  the  above  description 
of  the  fields.  Thus  the  bending  moment  Mr  and  the  shear  force  Sr  acting  on  the  tranverse  plane 
at  the  remote  location  can  be  accurately  measured  by  means  of  strain  gauges  on  the  surface  of  the 
specimen.  Denoting  the  value  of  the  moment  on  the  ligament  as  inferred  from  the  measurements 
at  the  remote  surface  location  by  Af£,  overall  quasi-static  equilibrium  requires  that 

Ml  =  Mr  +  SrLr  (3.6) 

where  Lr  is  the  distance  of  the  remote  section  along  the  beam  from  the  cracked  section.  Since  the 
remote  section  being  considered  is  about  the  mid-way  between  the  cracked  section  and  the  roller 
suPPort5  Mr  and  SrLr  will  be  of  comparable  magnitude. 

In  a  dynamic  fracture  experiment,  the  inertial  response  of  the  material  between  the  cracked 
section  and  the  remote  section  must  be  included  in  (3.6).  Such  inertial  contribution  has  been 
included  in  the  analysis  in  [1].  However  if  the  time  to  fracture  initiation  is  greater  than  the 
transition  time,  the  inertial  contribution  is  small  and  as  given  by  (3.6)  is  the  appropriate 
moment  for  use  in  the  computation  of  J^c. 
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4.  NUMERICAL  PROCEDURE 


The  full-held  3-D  numerical  simulation  has  several  objectives.  It  will  be  shown  that  under 
dynamic  conditions,  the  local  value  of  J  at  points  along  the  crack  front  can  be  accurately  calculated 
using  the  finite  domain  (volume)  integral  representation  (2.13).  We  then  examine  the  effect  of 
transient  loading,  geometry  and  plastic  deformation  on  the  variation  of  J  along  the  crack  front. 
The  usefulness  of  of  the  transition  time  concept  and  the  accuracy  of  the  deep-crack  J  formula 
(3.5)  based  on  measurable  load  and  deformation  histories  will  be  assessed.  Finally  the  suitablity  of 
proposed  surface  measurements  for  the  determination  of  J  will  be  clarified. 


4.1  Finite  element  procedure 

The  duration  of  the  event  of  interest  relative  to  the  time  for  wave  passage  through  the  structure 
is  the  important  consideration  in  the  choice  of  an  integration  scheme  for  a  dynamic  analysis.  Our 
primary  interest  is  the  response  of  a  cracked  specimen  under  impact  load  over  a  time  span  which 
is  large  compared  to  the  time  for  the  wave  to  travel  the  length  of  the  specimen.  However  we  are 
also  interested  in  the  dynamic  response  of  the  specimen  near  the  transition  time  when  the  discrete 
stress  waves  may  still  be  important.  An  accurate  resolution  of  fields  at  early  time  require  that  small 
time  steps  be  taken.  Comparisons  of  relative  efficiency  of  explicit  and  implicit  schemes  have  shown 
that  three-dimensional  numerical  simulations  of  dynamic  problems  can  be  more  econcomically 
carried  out  using  explicit  schemes  [18].  Even  if  the  higher  frequencies  are  not  of  interest,  the  large 
bandwidth  that  arises  from  a  3-D  model  and  the  overhead  connected  with  an  out-of-core  solver 
will  render  an  implicit  scheme  ineffective.  Guided  by  these  considerations,  we  employed  a  lumped 
mass  explicit  integration  scheme  for  all  our  calculations.  The  resulting  discrete  equations  of  motion 
are  uncoupled  and  the  nodal  accelerations  are  obtained  by  trivial  inversion  of  the  diagonal  mass 
matrix.  Our  numerical  experimentations  suggest  that  the  8-node  hexahedron  element  is  optimal 
for  the  intended  numerical  study.  All  numerical  results  reported  in  this  paper  were  obtained  using 
a  modified  version  of  the  finite  element  program  FEAP  [19]. 

To  alleviate  potential  numerical  difficulties  associated  with  deformation  into  the  fully  yielded 
plastic  state,  we  employ  the  so-called  B-bar  method  for  the  formation  of  the  stiffness  matrices 
[20,21],  In  the  case  of  the  8-node  trilinear  hexahedron  element,  the  volumetric  components  of  B 
are  obtained  by  1-point  quadrature  while  the  deviatoric  components  are  obtained  by  the  2x2x2 
quadrature  rule.  The  former  integration  results  in  a  uniform  hydrostatic  pressure  throughout  the 
element.  As  a  safequard  against  spurious  pressure  mode,  the  modified  B  matrix  [1]  is  the  normalized 
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sum  of  the  B  matrix  given  in  [21]  and  the  regular  B  matrix  [22]: 

B  =  Bdev  +  Bvo'  +  e(Bvol-Bvo1)  0<e<l  (4.1) 


The  regular  B  matrix  is  recovered  from  (4.1)  by  setting  e  —  1  while  the  B  matrix  in  [21]  is  obtained 
for  €  =  0.  In  our  analysis,  e  =  0.05  was  used. 


4.2  Discrete  formulas  for  the  evaluation  of  J 

The  discrete  form  of  J  (2.13)  based  on  2  X  2  X  2  Gaussian  quadrature  appropriate  to  the  8-node 
linear  hexahedron  element  is 


J 


£  t{ 

elements  p~  1 


all  elements  p= 
in  V 


dui  dqk  TTdqk 


aijdxk 


dxj 


U 


dxk 


+  P 


dui 


dt2  dxk 


Qk 


det 


dXr_ 

drjr, 


(4.2) 


Here  the  field  quantities  including  q and  its  spatial  derivatives  are  evaluated  at  quadrature  points 
(p  =  1,2,...,  8)  and  weighted  by  wp  and  the  determinant  of  the  Jacobian  matrix  det  (dxm/ dr]7n) 
[22].  The  acceleration  d2ui/dt 2  at  quadrature  points  are  obtained  from  nodal  values  using  the 
trilinear  interpolation  function.  To  maintain  consistency  with  the  isoparametric  formulation,  the 
values  of  q k  and  dq^/dxi  at  the  quadrature  points  are  evaluated  from 


8 

qk  =  ^  NiQki 

i=  i 


dqk 

dx{ 


8  3 


££ 


dNi  drjm 
drjm  dxi 


Qki 


(4.3) 


Here  Nj  is  the  trilinear  shape  function  and  Qfc/  is  the  value  of  qk  associated  with  the  Ith  local  node 
of  an  element.  Nodal  values  Qki  can  be  assigned  in  accordance  with  any  smooth  function  as  long 
as  the  appropriate  values  for  qk  on  St,  S0 ,  S+  and  S-  are  obtained  (see  Section  2.2).  Numerical 
experimentations  have  shown  that  the  value  of  J  is  insensitive  to  the  choice  of  Qki •  Thus  mesh 
design  and  convenience  are  the  only  considerations  in  selecting  smooth  functions  for  generating 
the  values  Qki-  A  more  detail  discussion  of  various  aspects  of  the  implementation  of  the  domain 
representation  of  J  has  been  given  in  [6,11]. 

The  value  of  J  at  a  nodal  point  (or  the  function  J(s)  along  the  crack  front)  can  be  determined 
in  the  following  manner.  We  consider  the  Mth  line  segment  which  includes  the  Mth  node  on  the 
crack  front.  The  line  segment  begins  and  ends  on  the  nodes  M  —  1  and  M  +  1  at  s*1  and  sj}1 
respectively.  Points  on  the  segment  advances  by  Am(s)  (in  the  normal  direction  in  the  plane  of  the 
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crack).  In  terms  of  the  present  variables,  the  energy  released  by  the  body  (2.7)  is  written  as 

sM 

Jm  -  fM  J(s)AM(s)ds  (4.4) 

''Sa 

Here  Jm  is  the  energy  release  computed  by  the  domain  integral  expression  (4.2).  We  assume  that 
the  local  energy  release  rate,  J(s),  is  a  continuous  function  of  s  and  can  be  represented  in  the  form 

NC 

j(s)  =  y,  N^s)  Jl  (4-5) 

L- 1 

where  Nc  is  the  total  number  of  nodes  on  the  entire  crack  front,  N l(s)  is  the  basis  function 
consistent  with  the  finite  element  shape  function  and  Jl  is  the  value  of  J  at  the  Lth  node  on 
the  crack  front.  Since  the  crack  front  region  is  modeled  by  8-node  trilinear  elements,  NL{s)  is  a 
piecewise  linear  function  of  s.  Similarly  the  function  A m(s)  has  the  form 

NC 

Am(s)  =  J]  Nl(s)  A &  (4.6) 

L-l 

where  Xfa  is  the  crack  advance  at  the  Lth  node  on  the  crack  front.  The  total  number  of  segments 
advanced  in  this  manner  is  equal  to  the  total  number  of  nodes  on  the  crack  front  Nc-  By  requiring 
be  equal  to  unity  for  L  =  M,  and  be  zero  for  L  ^  M,  (4.6)  simplifies  to 

A  m(s)  =  Nm(s)  (4.7) 


Substituting  (4.5)  and  (4.7)  in  (4.4),  we  have 


J  M 


Nl(s)Nm(s)  ds 


(4.8) 


We  repeat  the  above  procedure  for  every  node  on  the  crack  front  to  obtain  a  system  of  Nc 
linear  equations  for  Nc  unknown  Jl  ' s.  The  value  of  J l  is  given  by  the  solution  to  (4.8),  i.e., 

^  -  rs*f 

Jl  ~  ^2  ALMl  J  m  where  ALm  ~  f  m  Nl{s)Nm(s )  ds  (4.9) 

M=  1  ^ SCL 


Since  the  basis  function  Nl(s)  spans  over  three  nodes  on  the  crack  front,  Alm  has  a  bandwidth  of 
three.  With  values  of  Jl  in  hand,  the  value  of  J  at  any  point  s  on  the  crack  front,  is  obtained  from 
(4.5). 
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The  value  of  Jl  can  also  be  obtained  by  appealing  to  the  approximation  (2.10)  for  J(s).  With 
this  approximation,  the  J  at  a  node  is  uncoupled  from  the  J* s  at  the  other  nodes  on  the  crack 
front,  i.e., 

Nl(s)  ds  (4-10) 

L 
a 

If  J  m  varies  smootlily  along  the  crack  front,  Jl  as  computed  by  (4.10)  differs  negligibly  from 
its  more  accurate  value  determined  by  (4.9).  However  if  a  particular  nodal  value  is  evaluated 
inaccurately  (e.g.  the  J  value  at  a  node  on  a  free  surface),  the  inaccuracy  can  have  an  adverse 
effect  on  the  values  of  J  at  the  neighboring  nodes  if  (4.9)  is  used.  The  inaccuracy  is  confined  to 
the  particular  node  if  (4.10)  is  used. 

The  average  value,  pertaining  to  the  entire  crack  front,  is  obtained  by  integrating  J(s)  over 
the  entire  crack  front  and  dividing  the  resulting  value  by  the  total  length  of  crack  front.  This  value, 
denoted  by  /ave>  may  also  be  determined  directly  from  (4.2).  In  this  calculation,  a  unit  virtual 
extension  of  the  entire  crack  front  (take  A(s)  =  1  for  all  the  nodes  on  the  crack  front)  is  imposed 
and  the  resulting  value  of  J  is  divided  by  the  crack  front  length  to  give  Jave.  We  note  that  the  two 
methods  of  computing  «/ave  are  equivalent. 

4.3  Calculations  of  moments  and  rotations 

The  moment  carried  by  the  uncracked  ligament  Ml  is  required  in  the  deep-crack  formula  (3.5) 
for  estimating  J .  The  moment  is  calculated  from  the  forces  acting  on  the  plane  of  the  ligament. 
On  the  symmetry  plane  aq  =  0,  the  nodal  forces  (normal  to  the  symmetry  plane)  acting  on  nodes 
which  are  positioned  at  the  same  x\  coordinate  are  summed  (each  row  of  nodes  being  parallel  to 
the  a'3-axis)  to  give  the  resultant  normal  nodal  force  on  that  row.  This  is  repeated  for  every  row 
of  nodes  on  the  ligament  plane.  The  location  where  the  linearly  interpolated  resultant  nodal  force 
changes  its  sign  is  taken  as  the  effective  hinge  axis.  Suppose  the  effective  hinge  axis  (a  line  parallel 
to  the  x 3 -axis)  is  directed  along  x\  =  Iil .  The  moment  carried  by  the  ligament  plane,  0  <  x\  <  b 
and  —B/2  <  x$  <  B  j 2,  is  given  by 

nl 

Ml  =  2  f2'(x\  -  hL)  (4.11) 

i=l 

Here  Nl  is  the  total  number  of  nodes  on  the  lialf-ligament  plane,  f2%  is  the  nodal  force  and  aq* 
is  the  aq  coordinate  of  the  ith  ligament  node.  In  writing  (4.11)  we  have  taken  advantage  of  the 
symmetry  condition  on  the  plane  ;r3  =  0. 

The  inferred  moment  carried  by  the  ligament  (3.6)  is  determined  from  the  stresses  given  by 


Jl  =  J . 


L 
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the  strip  of  elements  near  the  free  surface  at  =  B/2  at  the  quarter-span.  These  element  values 
should  be  representative  of  surface  measurements  taken  in  experiments.  To  facilitate  the  moment 
computations,  the  elements  are  arranged  so  that  their  centroids  lie  on  the  transverse  plane  £2  —  Lr. 
The  moment  inferred  from  near  surface  stresses  is 


Mi 


r  . 

B  [^22lJixl  -  hR )  + 


3  =  1 


(4.12) 


where  B  is  the  thickness  of  the  specimen,  Nr  is  the  number  of  elements  in  the  surface  layer,  x(  is 
the  element  centroid’s  x\  coordinate  and  ZJ  is  the  length  of  the  element  in  the  aq-direction.  The 
normal  and  shear  stresses  <r/2  and  <tX2  at  the  element’s  centroid  are  obtained  by  averaging  the 
respective  values  at  the  eight  Gauss  points.  The  linearly  interpolated  normal  stress  changes  its  sign 
at  x\  =  Jir.  This  defines  the  hinge  axis  for  evaluating  the  Mr  component  (the  first  of  two  terms 
on  the  right  side  of  (4.12))  of  the  inferred  moment. 

The  rotation  of  the  cracked  section  about  the  hinge  axis,  Or,  is  given  by 

eL{t)  =  m /  (H  -  hL(t))  (4.13) 


Here  6  is  the  average  of  the  opening  displacements  of  the  nodes  at  the  crack  mouth  (along  xx  =  H ), 
and  Hr  is  the  hinge  position  (determined  in  connection  with  (4.11))  at  time  t .  It  may  be  noted 
that  the  opening  displacement  at  the  crack  mouth  vary  negligibly  with  £3.  An  inferred  rotation, 
which  is  very  useful  from  the  point  of  view  of  dynamic  fracture  experiments,  is  0*L  given  by 

n(t)  =  S(t)/  ( H-hl )  (4.14) 

Here  li*L  is  the  hinge  axis  given  by  the  slip-line  field  solution  in  [17].  In  particular  h*L/b  =  0.63 
which  is  very  close  to  the  hinge  axis  given  by  the  finite  element  solution  for  the  fully  yielded  state. 
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5.  FINITE  ELEMENT  ANALYSIS 


5.1  Computational  model 

The  relative  dimensions  of  the  three-point-bend  specimen  are  L/H  =  2.0,  a/H  =  0.5,  B/H  — 
0.5  and  Lt/H  =  2.25  where  H  is  the  width  of  the  specimen,  L  is  the  length  between  the  cracked 
section  and  the  remote  support,  a  is  the  crack  length,  B  is  the  thickness  and  Lt  is  the  half-length  of 
the  specimen  (see  Fig.  2a).  These  relative  dimensions  are  in  accord  with  ASTM  specifications  for 
the  specimen  for  quasi-static  fracture  toughness  evaluation.  A  typical  dynamic  fracture  specimen 
has  a  width  H  equal  to  2  inches  and  relative  crack  length  a/H  ranging  from  0.4  to  0.7  [2]. 

Only  a  quarter  of  the  three-point-bend  specimen  (region  x2  >  0,  £3  >  0)  need  be  modeled 
since  the  problem  possesses  reflective  symmetry  with  respect  to  the  crack  plane  (x2  =  0)  and 
mid-plane  (£3  =  0).  The  finite  element  model  of  the  quarter-specimen  constructed  with  8-node 
trilinear  hexahedron  (brick)  elements  is  shown  in  Fig.  2b.  Consider  a  layer  of  elements  (one  layer 
thick  in  the  £3  direction)  which  spans  the  x\-x2  plane.  The  element  size  in  the  £1-2*2  plane  is 
gradually  increased  with  element  distance  from  the  crack  plane.  The  identical  planar  mesh  is 
repeated  along  the  233-axis  to  form  six  layers  of  elements.  To  accomodate  the  anticipated  strong 
variations  of  stresses  and  deformation  (with  respect  to  £3)  as  the  free  surface  is  approached,  the 
thickness  of  the  elements  is  gradually  reduced.  The  six  layer  of  elements  are  bounded  by  the 
planes  £3 / B  —  0.0,  0.17,  0.27,  0.35,  0.41,  0.46,  0.5.  The  crack  front  is  surrounded  by  six  rings  of 
specially  arranged  crack-tip  elements.  Each  ring  consists  of  four  brick  elements  which  are  collapsed 
into  wedges.  The  outermost  ring  at  the  free  surface  is  shown  in  Fig.  2b.  Displacement  gradients 
within  these  elements  contain  terms  of  0(l/r)  in  the  £i-£2  plane  [23]  which  nearly  match  the 
near-tip  strain  singularities  associated  with  power-law  hardening  materials.  The  mesh  has  a  total 
of  1692  elements  (282  elements  on  each  layer)  and  2233  nodes. 

On  the  plane  of  the  uncracked  ligament  and  mid-plane,  normal  displacements  are  prescribed 
to  vanish.  The  model  is  roller-supported  along  the  line  X\  —  FT,  x2  ~  L.  A  uniform  pressure  is 
applied  over  a  quarter-region  bounded  by  0  <  £2  <  H/8  and  0  <  £3  <  J9/2  on  the  plane  xx  =  0 
so  that  the  quarter-load  is  Pf 4.  Twelve  element  surfaces  constitute  the  plane  of  contact.  At  very 
early  times,  corresponding  to  the  first  few  passes  of  the  discrete  waves  across  the  length  of  the 
specimen,  our  calculations  indicate  that  the  specimen  will  lift  off  from  its  roller  supports.  However 
the  imposed  roller  support  conditions  on  the  computational  model  did  not  allow  for  any  motion  in 
the  X\ -direction.  We  did  not  refine  our  computational  model  to  accomodate  this  phenomenon  since 
the  interval  of  interest  to  us  falls  far  beyond  the  time  where  discrete  wave  effects  are  important. 

The  material  behavior  is  taken  to  be  governed  by  the  J2  flow  (rate  independent)  theory  of 
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plasticity.  Under  uniaxial  stressing  the  material  deforms  according  to 


e/eo  =  0-/uo,  e<e0 
e/fo  =  (o-/o-o)n,  e  > 


(5.1) 


where  cr0  and  e0  are  the  yield  stress  and  strain  related  by  a„  =  Ee0,  E  is  Young’s  modulus  and  n 
is  the  strain  hardening  exponent.  In  our  calculations,  we  have  taken  n  to  be  10  and  Poisson’s  ratio 
v  to  be  0.3,  these  values  being  representative  of  high-strength  structural  steels. 


5.2  Dynamic  Analysis 

To  obtain  an  estimate  of  the  stable  time  step  for  the  mesh  depicted  in  Fig.  2b,  we  obtained  the 
largest  eigenvalue  of  the  global  elastic  stiffness  matrix.  From  this  value,  the  dimensionless  critical 
time  step  was  determined  to  be  A t  c\/ H  —  0.014.  Uniform  time  increments  of  At  c\/ H  —  0.012  was 
chosen  for  the  entire  simulation.  The  load  applied  at  the  mid-span  increases  linearly  in  the  interval 
0  <  tci/H  <  40,  and  thereafter  is  kept  constant  at  PL/2Mo  =  1.2  (see  insert  in  Fig. 6).  This 
loading  rate  is  representative  of  the  nearly  constant  load  application  rate  that  has  been  achieved 
in  three-point-bend  bar  fracture  experiments  through  proper  choice  of  impact  absorber  material 
[2].  Altogether  12000  uniform  time  increments  were  taken  and  the  value  of  J  and  pertinent  field 
quantities  were  evaluated  at  every  50  increments  corresponding  to  a  time  interval  of  tc\/H  —  0.6. 
For  this  particular  simulation,  the  CPU  time  expended  on  the  Cray-XMP  is  about  six  hours. 

We  extracted  /ave  and  the  value  of  J  at  each  crack  front  node  from  finite  element  fields  using 
the  procedure  discussed  in  Section  4.2.  In  the  J  calculations,  values  of  qi  and  q 3  are  prescribed  to  be 
zero  everywhere  in  the  domain  since  the  normal  to  the  (straight)  crack  front  is  in  the  x\ -direction. 
To  calculate  J  m  using  (4.2),  q\  is  set  equal  to  unity  at  the  Mth  node  and  zero  at  all  other  nodes 
on  the  crack  front.  Within  a  domain  chosen  for  each  calculation,  the  values  of  q\  are  in  accord 
with  the  conditions  given  in  Section  2.2  and  further  elaborated  upon  in  Section  4.2.  Four  different 
domains  were  employed  for  the  evaluation  of  J  m-  The  smallest  domain  had  twenty  elements  while 
the  largest  domain  had  more  than  four  hundred  elements.  Typically  the  value  of  J  m  associated 
with  a  domain  is  within  two  percent  of  its  mean  value.  The  near  domain  independence  of  J  is  an 
indication  of  the  quality  of  the  numerical  solution.  The  mean  value  of  J  is  used  in  (4.9)  and  (4.10) 
for  the  evaluation  of  J  at  the  nodes  and  J(s)  at  the  point  s  on  the  crack  front.  We  computed 
Jave  by  integrating  J(s)  along  the  crack  front  and  then  divide  the  result  by  the  crack  length.  As 
a  check,  the  average  value  is  computed  using  (4.2)  by  taking  qi  to  equal  unity  on  all  crack  front 
nodes.  The  Jave  values  determined  by  the  two  methods  are  exactly  equal  as  they  must  be. 
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The  increase  of  JwJ<r0e0a  with  time  tci/H  is  plotted  in  Fig.  4.  Also  shown  in  the  figure, 
is  J  at  the  symmetry  plane  (z3 /B  =  0)  and  at  the  free  surface  (x^/B  —  0.5),  denoted  by  Jm\<\ 
and  Jedge  respectively.  The  local  values  at  these  two  locations  and  all  other  crack  front  locations 
are  determined  using  (4.10).  Shown  in  Fig.  5  is  the  variation  of  J  along  the  crack  front  at  four 
different  times.  The  local  J  values  at  each  time  have  been  normalized  by  the  corresponding  JaVe* 
The  increasing  through-thickness  variation  of  J  with  the  progress  of  time  reveals  the  effect  of 
plasticity  which  dominates  the  long  time  response  of  the  specimen.  Similar  plasticity  effects  on 
three  dimensional  near  tip  fields  have  been  observed  in  a  quasi- statically  loaded  specimen  [4], 

The  actual  moment  carried  by  the  ligament  Ml  is  calculated  by  (4.11).  For  subsequent  dis¬ 
cussion,  we  normalize  Ml  by  the  limit  moment  given  by  the  plane  strain  rigid-perfectly  plastic 
analysis  [17].  The  limit  moment,  M0  has  the  value  0.364  cr0b2B  where  aQ  is  the  yield  stress,  b  is 
the  ligament  length  and  B  is  the  thickness  of  specimen.  The  computed  location  of  the  effective 
hinge  axis  in  the  fully  yielded  state  is  very  close  to  the  location  given  in  [17].  Plotted  in  Fig.  6 
is  the  increase  of  the  normalized  moment  ML/M0  with  dimensionless  time  tci/H .  The  effect  of 
inertial  screening  is  clearly  evident  at  tc\/H  <  20.  For  times  tc\/H  greater  than  about  20,  the 
moment  carried  by  the  ligament  rises  rapidly.  Beyond  tc\/H  of  about  40,  Ml  begins  to  level  off 
and  slowly  approach  the  equilibrium  moment.  The  rotation  of  the  cracked  section  9l  with  time  is 
nearly  identical  to  that  observed  in  the  2-D  simulation  of  the  same  impact  problem  [1]  and  will  not 
be  discussed. 

Through- thickness  variation  of  the  near  tip  stresses  provide  further  insight  into  the  3-D  nature 
of  the  fields.  We  examine  the  stresses  on  two  lines  parallel  to  the  &3-axis  and  located  near  to  the 
crack  front.  Of  particular  interest  in  fracture  initiation  studies  is  the  tensile  stress  (722*  These 
stresses  are  evaluated  at  element  centroids  located  at  r/H  =  0.035,  9  =  27.3°  and  r/H  —  0.099,  9  — 
18.4°  where  r  and  9  are  the  polar  coordinates  centered  at  the  crack  front  and  9  is  measured  from 
the  crack  plane.  The  through-thickness  variation  of  (T22  (normalized  by  the  yield  stress  aG)  at  four 
different  times  along  the  two  line  locations  is  shown  in  Fig.  7.  At  shorter  times  when  the  material 
responds  elastically,  the  fields  may  be  characterized  as  2- dimensional.  At  longer  times  dominated 
by  the  plastic  response  of  the  material,  the  stresses  develop  strong  through-thickness  variation. 
This  trend  is  similar  to  the  variation  of  J  along  the  crack  front  (see  Fig.  5). 

At  the  remote  plane  located  approximately  mid- way  between  the  crack  plane  and  roller  sup¬ 
ports  ( X2  =  L/2)  there  is  hardly  any  variation  of  the  axial  stresses  and  the  shear  stresses  through 
the  thickness  at  all  times,  i.e.,  the  fields  may  be  characterized  as  plane  stress.  Furthermore  the 
distribution  of  the  normal  strain  and  stress  with  .xi  is  very  nearly  in  agreement  with  the  distribu¬ 
tion  according  to  elementary  beam  theory.  It  follows  that  the  load  carried  by  the  ligament  can  be 
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inferred  from  surface  measurements  taken  in  the  vicinity  of  the  quarter-span.  This  is  significant 
since  this  load  is  the  input  measurements  for  the  application  of  the  deep-crack  formula  for  J .  This 
will  be  taken  up  in  the  next  section. 

In  the  interval  0  <  t  <  tr,  the  material  response  is  elastic.  At  the  transition  time  tT  elements 
immediately  adjacent  to  the  crack  tip  are  partially  yielded.  As  discussed  in  Section  3.2  it  is  doubtful 
that  at  this  stage  of  the  transient  the  near-tip  fields  are  sufficiently  established  so  that  they  can  be 
characterized  by  J.  At  larger  times  the  plastic  zone  emanating  from  the  crack  tip  spreads  across 
the  ligament  and  joins  up  with  the  plastic  zone  emanating  from  the  plane  of  impact.  Effective 
stress  contours  at  tc\/H  =  144  are  shown  in  Fig.  85  the  contour  labeled  by  1.0  is  the  boundary 
of  the  plastic  zone.  Plastic  zone  shape  in  the  x^x2  plane  resembles  the  quasi-static  plane  strain 
plastic  zone  [16,23].  The  extent  of  the  plastic  zone  near  the  free  surface  ( xz/B  =  ±0.5)  is  larger 
than  that  at  the  mid-plane  {xz/B  =  0)  of  the  specimen. 

To  gain  further  insight  into  the  behavior  of  the  dynamically  loaded  specimen,  results  from 
the  dynamic  analysis  are  presented  in  a  way  that  will  allow  a  direct  comparison  with  the  quasi¬ 
static  response  of  the  specimen.  First  a  quasi-static  analysis  based  on  the  identical  mesh  is  briefly 
discussed.  A  small  load  increment  is  applied  so  that  the  Gauss  points  nearest  to  the  crack  tip  are 
on  the  verge  of  yielding.  Thereafter  the  load  increment  is  adjusted  so  that  the  plastic  zone  increases 
gradually.  A  total  of  70  increments  were  taken  to  reach  the  final  load  of  Ml/Mq  =  1.3.  On  the 
average  seven  Newton- Raphson  iterations  were  required  to  bring  the  imbalanced  residual  forces  to 
a  normalized  error  level  of  10-5. 

Values  of  dL /e0  and  Jave/a0e0a  as  functions  of  ML/M0  as  determined  from  the  dynamic  analysis 
are  shown  as  solid  line  curves  in  Fig.  9a  and  9b.  The  triangles  are  the  results  from  a  quasi-static 
analysis.  The  close  agreement  is  indeed  remarkable.  Though  not  detectable  in  the  plots,  the 
quasi-static  and  dynamic  values  differ  by  about  20  percent  at  load  level  of  Ml/M0  —  0.1;  at 
even  earlier  times  and  lower  loads  the  differences  are  larger.  Plots  of  J / cr060a  versus  0lI6o  f°r 
both  dynamic  and  quasi-static  cases  are  shown  in  Fig.  9c.  We  observe  that  with  the  appropriate 
normalizations,  the  relationship  between  normalized  Jave  and  normalized  displacement  (or  rotation 
of  the  ligament),  appear  to  be  insenitive  to  loading  rate  in  the  regime  dominated  by  structural 
inertial  and  deformation  energy.  A  similar  observation  can  be  made  on  the  relationship  between  J 
and  the  moment. 

5.3  Comparison  of  (precise)  J  and  Jdc 

With  the  moment  Ml  (4.11)  and  the  rotation  0L  (4.13)  as  the  input  measurements  in  (3.5), 
the  value  of  Jdc  as  a  function  of  time  is  obtained.  For  the  purpose  of  appraising  the  accuracy  of  the 
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deep-crack  J  formula,  the  value  of  Jdc  is  divided  by  the  average  value,  Jave.  The  latter  is  the  precise 
value  determined  from  the  computed  field  quantities  using  the  domain  integral  expression  (4.2).  A 
ratio  of  unity  represents  perfect  agreement  between  the  deep-crack  estimate  and  the  precise  value 
Jave.  The  ratio  Jdc/J*ve  is  plotted  against  tci/H  in  Fig.  10a.  During  the  early  transient,  the 
ratio  is  substantially  larger  than  unity,  then  falls  below  unity  and  starts  to  increase  at  about  the 
transition  time.  It  levels  off  at  a  value  of  1.06  at  about  twice  the  transition  time. 

The  ratio  as  a  function  of  the  actual  moment  carried  by  the  ligament  is  shown  by  the  solid 
line  curve  in  Fig.  10b.  We  have  noted  that  the  actual  moment  carried  by  the  ligament  cannot 
be  measured  in  an  experiment,  and  a  measurable  moment  was  introduced  Section  3.3.  For  this 
particular  calculation,  the  measurable  moment  is  determined  by  the  fields  at  the  free  surface  near 
to  the  quarter-span,  x2  =  1.068#.  As  discussed  earlier  the  moment  and  shear  force  can  be  reliably 
determined  from  stresses  near  the  surface  at  the  chosen  location  since  through-thickness  variation  of 
the  fields  near  the  quarter-span  is  negligible  and  the  aq -distribution  of  the  numerically  determined 
fields  is  nearly  exactly  given  by  elementary  beam  theory.  Using  the  inferred  moment  from  (3.6) 
and  Q*l  from  (4.14)  in  the  deep-crack  formula  (3.5),  a  ‘measurable’  Jdc  is  obtained  .  This  estimate 
of  J  normalized  by  the  precise  /ave  is  shown  by  the  dashed  line  curve  in  Fig.  10b.  The  agreement 
between  the  dashed  line  curve  based  on  the  ‘measurable’  moment  and  rotation,  and  the  solid  line 
curve  based  on  the  actual  moment  and  rotation  is  quite  good. 
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6.  DISCUSSION 


It  has  become  fairly  common  practice  to  write  the  relationship  between  J  and  the  work  done 
on  a  specimen  by  the  applied  loads  W  as 

Javc  =r]Wl{Bb)  (6.1) 

where  b  is  the  length  of  uncracked  ligament  and  Bb  is  the  area  of  the  uncracked  ligament.  In 
writing  (6.1),  we  have  followed  the  notation  in  this  paper;  in  the  published  literature  (particularly 
papers  on  fracture  toughness  measurements,  e.g.  [2])  the  left  hand  side  of  (6.1)  is  usually  written 
as  J.  Relationship  (6.1)  is  a  quasi-static  loading  deformation  theory  result  and  should  provide  an 
accurate  estimate  of  J  if  proportional  loading  is  nearly  satisfied.  For  a  deeply  cracked  bar,  the 
values  of  ‘eta  factor’  are  available  for  remote  loading  ranging  from  pure  stretch  to  pure  bending 
[24].  The  significance  of  the  ‘eta  factor’  is  that  it  permits  the  the  determination  of  J  directly  from 
the  work  done  on  the  specimen  by  the  applied  load.  Guided  by  (6.1),  we  introduce  a  similar  ‘eta 
factor’  for  transient  loading  defined  by 

■n  =  J^Bb/  J ^  ML(t')dOL(t')  (6.2) 

where  the  integral  in  (6.2)  is  the  work  done  on  the  ligament  by  the  load  Ml  carried  by  the  ligament 
(see  Section  3.2).  The  variation  of  rj  with  normalized  moment  is  shown  in  Fig.  11  for  times  greater 
than  the  transition  time.  The  solid  line  curve  is  determined  from  the  results  for  present  dynamic 
problem  using  (6.2).  The  triangles  are  values  determined  from  finite  element  results  for  the  quais- 
statically  loaded  specimen.  While  ‘eta  factor’  varies  negligibly  with  the  applied  moment  in  the 
quasi-static  case  (as  it  must),  the  value  for  the  dynamic  case  starts  at  about  2.17  and  decreases  to 
about  1.95  for  Ml/Mq  equal  to  1.2. 

The  transition  time  concept  was  introduced  in  [1]  to  provide  an  estimate  of  the  range  of  time 
during  which  the  crack  tip  field  stabilizes  and  the  deep-crack  formula  is  applicable  under  transient 
loadings.  In  this  paper,  we  have  shown  that  the  estimate  of  the  transition  time  based  on  an 
elementary  2-D  model  analysis  is  in  good  agreement  with  the  transition  time  obtained  from  the 
full-field  3-D  transient  solution  (see  Fig.  3).  The  usefulness  of  the  concept  for  the  particular  class 
of  dynamic  fracture  experiments  being  discussed  is  apparent  from  the  plots  in  Figs.  10  and  11. 

Full-field  finite  element  simulation  of  a  standard  ASTM  three-point-bend  bar  with  relative 
crack  depth  0.5  yields  a  transition  time  of  about  24 Hjc0  (or  28 H/ C\).  This  time  is  roughly  equal 
to  the  time  required  for  the  longitudinal  wave  to  make  six  traverses  of  the  span  of  the  specimen, 
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or  for  the  slower  shear  wave  to  make  three  passes.  Indeed  full-field  calculations  indicate  that  the 
fields  have  stabilized  at  times  larger  than  twice  the  transition  time.  At  such  times  the  total  kinetic 
energy  is  less  than  20  percent  of  the  total  stress  work  on  the  body  (see  Fig.  3). 

Similar  observations  concerning  the  transition  from  a  response  dominated  by  individual  stress 
waves  to  a  response  dominated  by  the  fundamental  structural  mode  has  also  been  discussed  in  a 
review  paper  by  Ireland  [25].  Ireland  introduced  an  effective  specimen  inertial  oscillation  period, 
r,  and  cited  numerous  experimental  data  which  showed  that  inertial  effects  are  dominant  for  times 
smaller  than  2  r.  His  proposed  empirical  relation  for  r  can  be  arranged  in  the  form 

r  =  1.68  x  \/2  S  —  (6.3) 

cQ 

where  S  is  defined  by  (3.2).  Based  on  the  relative  specimen  dimensions  employed  in  the  present 
analysis,  a  value  of  r  equal  to  23.8 H / cQ  is  obtained  from  (6.3)  which  is  remarkably  close  to  the 
value  we  got  from  the  finite  element  solution  and  the  model  analysis  (3.1).  Kalthoff  [26]  has  also 
noted  that  elastic  near-tip  fields  at  short-times  cannot  be  predicted  by  a  static  analysis.  In  the 
short-time  transient  regime  (t  <  3  r),  Kalthoff  and  coworkers  (see  publications  referenced  in  [26]) 
advocate  the  use  of  impact  response  curves  for  estimating  the  dynamic  stress  intensity  factor. 

This  study  is  directed  at  ductile  materials  where  fracture  initiation  takes  place  after  extensive 
yielding  of  the  ligament  has  occured.  Our  analysis  suggests  that  interpretable  fracture  toughness 
data  can  be  obtained  from  dynamically  loaded  three-point-bend  steel  specimen  if  the  time  to 
fracture  initiation  is  larger  than  about  twice  the  transition  time.  In  particular,  the  transition  time 
for  an  ASTM  three-point-bend  steel  specimen  of  width  5  cm  (this  dimension  being  typically  denoted 
by  W  in  the  standards)  and  crack  length  to  width  ratio  between  0.5  to  0.7  ranges  from  about  250  to 
400  (is.  Thus  the  time  to  fracture  initiation  should  be  greater  than  about  500  to  800  (is.  Fracture 
initiation  in  ductile  materials  at  much  shorter  times  (less  than  30  ( is )  can  be  achieved  in  the  stress 
wave  loaded  cracked  round  bar  experiment  of  Costin,  Duffy  and  Freund  [27].  A  full-field  transient 
finite  element  analysis  of  the  fracture  experiment  has  been  carried  out  by  Nakamura,  Shih  and 
Freund  [28]. 
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FIGURE  CAPTIONS 


Fig.  1  (a)  Crack  tip  contour  T  inscribed  on  the  plane  locally  perpendicular  to  the  crack  front  where 
s  is  the  arclength,  (b)  virtual  crack  advance  within  sa  and  sb  represented  by  Z*.(s),  (c)  unit 
normal  nif  pointing  ‘away’  from  volume  enclosed  by  <S+,  S-  and  SQ . 

Fig.  2  (a)  Schematic  of  three-point-bend  specimen.  The  global  coordinates  Xi ,  x2  and  are  indi¬ 
cated.  (b)  Finite  element  model  of  quarter- specimen. 

Fig.  3  Ratio  of  total  kinetic  to  deformation  energies  as  determined  from  an  elementary  elastic  model 
analysis.  For  comparison  the  finite  element  results  are  included.  The  transition  time  t t  is 
indicated. 

Fig.  4  Behavior  of  Jave,  and  local  J1 s  at  the  mid-point  and  the  edge  of  the  crack  front. 

Fig.  5  Variation  of  J  along  crack  front  at  four  different  times.  The  pointwise  J  values  are  normalized 
by  the  corresponding  Jave. 

Fig.  6  Moment  carried  by  ligament  plane  as  a  function  of  elapsed  time.  The  applied  load  is  shown  in 
insert. 

Fig.  7  Through-thickness  variations  of  normal  tensile  stress  (722*  The  locations  of  the  two  lines  (par¬ 
allel  to  the  £3-axis  are  (a)  x\/H  =  0.469,  x2 /H  =  0.016  ( r/II  —  0.035,  9  —  27.3°),  (b) 
Xi/H  =  0.406,  x2/H  =  0.031  (r/H  =  0.099,  9  =  18.4°). 

Fig.  8  Effective  stress  contours  for  cre/(70  =  0.8,  1.0,  1.2,  1.3,  at  tcx/H  =  144. 

Fig.  9  Relationships  between  moment  carried  by  ligament,  rotation  of  cracked  section  and  JaVe-  Static 
results  shown  by  the  triangles  are  included  for  comparison. 

Fig.  10  Comparison  of  J^c  with  J&ve  (a)  as  a  function  of  time,  (b)  as  a  function  of  moment  carried  by 
ligament.  The  dashed  line  is  based  on  measurable  moment  and  rotation. 

Fig.  11  Dependence  of ‘eta  factor’  on  moment  carried  by  ligament  for  times  greater  than  the  transition 
time.  Values  for  the  quasi-statically  loaded  specimen  are  shown  by  the  triangles. 
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